#include "mathutil.h"
#include "lm.h"

using namespace std;

float LanczosSinc(float t, float tau ) {
	t = std::abs(t);
	if (t < 1e-6)
		return 1.0f;
	else if (t > 1.0f)
		return 0.0f;
	t *= M_PI;
	float sincTerm = std::sin(t*tau)/(t*tau);
	float windowTerm = std::sin(t)/t;
	return sincTerm * windowTerm;
}

bool SolveLinearSystem2x2(const float A[2][2],
	const float B[2], float *x0, float *x1) {
		float det = A[0][0]*A[1][1] - A[0][1]*A[1][0];
		if (fabsf(det) < 1e-10f)
			return false;
		*x0 = (A[1][1]*B[0] - A[0][1]*B[1]) / det;
		*x1 = (A[0][0]*B[1] - A[1][0]*B[0]) / det;
		if (isnan(*x0) || isnan(*x1))
			return false;
		return true;
}

bool Quadratic(float A, float B, float C, float *t0, float *t1) {
	float discrim = B * B - 4.f * A * C;
	if(discrim < 0.)
		return false;
	float rootDiscrim = sqrtf(discrim);
	float q;
	if (B < 0) q = -0.5f * (B - rootDiscrim);
	else q = -0.5f * (B + rootDiscrim);
	*t0 = q / A;
	*t1 = C / q;
	if(*t0 > *t1)
		swap(*t0, *t1);
	return true;
}

float Polar(float x, float y) {
	float d = atanf(x/y);
	if(d < EPSILON(6))
		d = 0;
	if(x > 0.) {
		if(y < 0.)
			d += M_PI;
	}
	else {
		if(y > 0.)
			d += 2 * M_PI;
		else
			d += M_PI;
	}
	return d;
}